#include<../data_tools/usefull_routines.c>

int main(int argc, char *argv[])
{
  int i,j,npart,nbines,ndisk,limpart;
  double *radio,rmin,rmax,deltar,rini,rfin;
  double vR,vang,vz,meanvR,meanvang,meanvz;
  double disvR,disvang,disvz;
  size_t *ind;
  char *infile,outfiles[1000];
  particulas *diskpart;
  infile = argv[1];
  nbines = atoi(argv[2]);
  limpart = 100;

  FILE *diskfile;
  diskfile = fopen(infile,"r");
  if(diskfile==NULL)
    {
      printf("No existe el archivo %s\n",infile);
    }
  ndisk = contadorParticulas(infile);

  diskpart = (particulas *)malloc((size_t)ndisk*sizeof(particulas));
  if(diskpart == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  radio = (double *)malloc((size_t)ndisk*sizeof(double));
  if(radio == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  ind = (size_t *)malloc((size_t)ndisk*sizeof(size_t));
  if(ind == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=0; i<ndisk; i++)
    {
      fscanf(diskfile,"%d %d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e",
	     &diskpart[i].id,
	     &diskpart[i].type,
	     &diskpart[i].pos[0],&diskpart[i].pos[1],&diskpart[i].pos[2],
	     &diskpart[i].vel[0],&diskpart[i].vel[1],&diskpart[i].vel[2],
	     &diskpart[i].mass,
	     &diskpart[i].pot,
	     &diskpart[i].posc[0],&diskpart[i].posc[1],&diskpart[i].posc[2],
	     &diskpart[i].velc[0],&diskpart[i].velc[1],&diskpart[i].velc[2],
	     &diskpart[i].mang[0],&diskpart[i].mang[1],&diskpart[i].mang[2],
	     &diskpart[i].mr,
	     &diskpart[i].r,
	     &diskpart[i].vc,
	     &diskpart[i].epsilon,
	     &diskpart[i].Ebinding,
	     &diskpart[i].je[0],&diskpart[i].je[1],&diskpart[i].je[2]);
      radio[i] = diskpart[i].posc[0];
      ind[i] = i;
      //printf("%e\n",diskpart[i].pot);
    }

  gsl_sort_index(ind,radio,1,ndisk);
  
  rmin = log10(radio[ind[0]]);
  //rmax = log10(radio[ind[ndisk-1]]);
  rmax = log10(10.0);
  printf("Done.\n");
  printf("\n");

  //Computing radial density profile
  printf("Computing dispertions velocities\n");
  FILE *dispertionfile;
  sprintf(outfiles,"%s_dispertion.dat",infile);
  dispertionfile = fopen(outfiles,"w");

  meanvR =  meanvang =  meanvz = 0.0;
 
  deltar = (rmax-rmin)/nbines;
  rini = rmin;
  for(i=0; i<nbines; i++)
    {
      rfin = rini+deltar;
      npart = 0;
      vR = vang = vz = 0.0;
      disvR = disvang = disvz = 0.0;

      for(j=0; j<ndisk; j++)
	{
	  if((radio[ind[j]]>=pow(10,rini)) && (radio[ind[j]]<pow(10,rfin)))
	    {
	      // if((diskpart[ind[j]].pos[2]>=-0.5) && (diskpart[ind[j]].pos[2]<0.5))
	      //	{
	      npart = npart+1;
	      vR = vR+diskpart[ind[j]].velc[0];
	      vang = vang+diskpart[ind[j]].velc[1];
	      vz = vz+diskpart[ind[j]].velc[2];
	      //	}
	    }
	}
      if(npart>=limpart)
	{ 
	  meanvR = vR/npart;
	  meanvang = vang/npart;
	  meanvz = vz/npart;
	}
      npart = 0;
      for(j=0; j<ndisk; j++)
	{
	  if((radio[ind[j]]>=pow(10,rini)) && (radio[ind[j]]<pow(10,rfin)))
	    {
	      //if((diskpart[ind[j]].pos[2]>=-0.5) && (diskpart[ind[j]].pos[2]<0.5))
	      //	{
	      npart = npart+1;
	      
	      disvR = disvR + pow((diskpart[ind[j]].velc[0]-meanvR),2);
	      disvang = disvang + pow((diskpart[ind[j]].velc[1]-meanvang),2);
	      disvz = disvz + pow((diskpart[ind[j]].velc[2]-meanvz),2);
		  //	}
	    }
	}
      if(npart>=limpart)
	{ 
	  disvR = disvR/(npart-1);
	  disvang = disvang/(npart-1);
	  disvz = disvz/(npart-1);
	  fprintf(dispertionfile,"%e %e %e %e %e\n",pow(10,rini),sqrt(disvR),sqrt(disvang),sqrt(disvz),meanvz);
	}
      rini = rfin;
    }
  printf("Done.\n");
  printf("\n");

  free(diskpart);
  free(radio);
  free(ind);
  fclose(diskfile);
  fclose(dispertionfile);

  return 0;
}
